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1.  Introduction 

Two  separate  themes  from  the  Computer  Vision  literature  come  together  in  this  paper:  the  use  of  rota- 
tinnally  symmetric  operators,  and  the  idea  that  several  modules  of  visual  perception  require  that  the  "most 
conservative"  solution  that  meets  a  given  set  of  boundary  conditions  be  computed.  The  two  themes  are 
combined  in  an  investigation  of  which  operator  to  use  in  the  interpolation  of  smooth  surfaces  from  one¬ 
dimensional  boundary  constraints.  Such  constraints  arise  naturally  in  a  variety  of  visual  problems. 

In  the  next  section  we  review  the  role  of  rotationally  symmetric  operators  in  Computer  Vision,  and  we 
derive  conditions  which  linear  and  quadratic  forms  in  the  first  and  second  directional  derivatives  must  satisfy 
in  order  to  be  rotationally  symmetric.  We  then  discuss  the  idea  that  vision  is  a  conservative  process,  citing 
examples  from  both  figure  perception  and  scene  analysis.  The  "most  conservative"  solution  is  modelled  using 
the  calculus  of  variations  to  find  the  minimum  function  that  satisfies  a  given  performance  index.  A  major 
problem  associated  with  the  use  of  the  calculus  of  variations  is  guaranteeing  the  existence  of  an  minimum 
function  (see  for  example  Courant  and  Hilbert  1953,  p.173).  A  theorem  of  Grimson(1981,  theorem  2)  proves 
that  a  sufficient  condition  for  the  existence  of  a  minimum  is  that  the  performance  index  should  be  a  seminorm 
on  the  space  of  functions.  The  condition  is  not  necessary.  For  example,  Horn(]981)  has  determined  the  curve 
that  minimizes  ti  e  in  'cgral  square  curvature  subject  to  tangency  conditions  at  the  end  points',  the  performance 
index  is  not  a  seminorm. 

Grimson(1981)  notes  that  many  intuitively  plausible  performance  indices  based  on  mean  and  Gaussian 
curvature  are  not  seminorms,  but  that  the  square  Tapiacian  f\x  +  2fJyy  +  f2yy  and  the  quadratic  variation 
fix  +  2  fly  -f  f2yv  arc.  Wc  show  here  that  any  quadratic  form  in  /„,  fiy ,  and  fyy  is  a  seminorm. 

To  further  constrain  the  choice  of  performance  index  in  the  infinite  set  of  quadratic  forms,  wc  require 
in  addition  that  the  quadratic  form  should  be  rotationally  symmetric.  We  prove  that  there  are  essentially  two 
different  choices:  the  square  Tapiacian  and  the  quadratic  variation.  All  the  remaining  possibilities  arc  linear 
combinations,  that  is,  form  a  vector  space  with  these  two  as  a  basis, 

To  choose  between  the  square  Tapiacian  and  the  quadratic  variation,  wc  consider  their  respective  Kuler 
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conditions  and  natural  boundary  conditions  (Courant  and  Hilbert,  1953).  The  Euler  conditions  are  identical, 
but  the  natural  boundary  conditions,  which  are  derived  from  the  statics  of  a  deformed  thin  plate,  favor  the 
quadratic  variation  since  they  offer  tighter  constraint  in  this  case. 

2.  Rotationally  symmetric  operators  in  vision 

A  major  concern  of  Computer  Vision  is  the  isolation  of  constraints  that  combine  with  the  information 
provided  in  the  image  to  yield  an  interpretation.  Early  work  on  polyhedra  (Clowes  1971,  Huffman  1971, 
Mackworth  1973,  Waltz  1972,  Sugihara  1978, 1981,  Kanadc  1981)  focussed  upon  the  discovery  of  constraints 
deriving  from  the  image  forming  process,  constraints  that  relate  image  fragments,  like  junctions  and  lines,  to 
their  scene  counterparts,  vertices  and  edges.  As  Computer  Vision  turned  its  attention  away  from  plane-faced 
objects  to  the  natural  world,  other  constraints  were  required.  Often  the  constraints  expressed  some  facet  of  the 
intuitive  notion  of  "smoothness"  and  did  so  in  a  way  that  supported  useful  computations  (Strat  1979,  Brooks 
1979.  Ikcuchi  and  Horn  1981,  Woodham  1978,  Horn  and  Schunck  1981).  Recently,  smoothness  and  image 
forming  have  been  combined  using  differential  geometry  (Grimson  1981,  Witkin  1981,  Binford  1981). 

One  constraint  that  is  usually  implicit,  but  is  occasionally  made  explicit,  expresses  the  idea  that  perceptual 
processes  arc  often  approximately  isotropic.  It  seems  that  humans  usually  do  not  show  strong  directional 
preferences  when  detecting  edges,  motion,  or  reflectance  boundaries.  We  seem  to  be  equally  adept  at  per¬ 
ceiving  the  layout  and  orientation  of  a  visible  surface  regardless  of  its  orientation  relative  to  the  view  vector. 
Ullman(1976)  ar3"?s  for  an  explicit  isotropy  constraint  in  his  work  on  subjective  contours  (sec  also  Knuth 
1979). 

Processes  that  arc  isotropic  arc  naturally  computed  by  rotationally  symmetric  operators,  since  the  values 
they  return  are  unaffected  by  the  coordinate  system  chosen  for  the  image.  Conversely,  rotationally  symmetric 
operators  compute  isotropic  information.  As  we  shall  see.  many  operators  that  have  been  proposed  for  vision 
arc  not  rotationally  symmetric  but  directionally  selective.  Some  authors  have,  however,  proposed  rotationally 
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Precise  definitions  of  rotational  symmetry  for  functions,  operators  (or  functionals),  and,  by  specialization, 
matrices  are  given  in  the  following  section.  In  the  rest  of  this  section  we  assume  that  'he  definitions  are  already 
understood. 

Some  kinds  of  blurring  in  an  image  forming  system  can  be  approximated  by  convolution  with  a 
Gaussian.  The  rotationally  symmetric  Gaussian  can  be  defined  by: 

G(r)  - 

Pratt(1978)  presents  several  techniques,  such  as  convolution  with  the  generalized  inverse  of  the  blur 
function,  for  restoring  the  image,  (see  for  example,  his  figures  14.2.1, 14.3.2). 

The  Laplacian  A  —  fvy  is  well  known  to  be  rotationally  symmetric1  and  its  use  has  been  proposed 
several  times  in  Computer  Vision  and  Image  Processing.  If  an  image  is  blurred  in  a  way  that  can  be  ap¬ 
proximately  modelled  by  passing  the  image  through  a  system  with  a  Gaussian  point  spread  function,  then  it 
can  be  sharpened  by  subtracting  a  multiple  of  its  Laplacian  (Rosenfcld  and  Kak  1976,  p.184,  Prewitt  1970,  p. 
107).  Pratt(1978,  figure  17.4.5)  illustrates  the  use  of  the  laplacian  for  enhancing  the  edges  in  an  image. 

Wcska,  Dyer  and  Roscnfcld(1976)  note  that  convolving  a  step  edge  with  a  Laplacian  operator  gives  rise  to 
a  pulse  pair:  a  negative  pulse  at  the  transition  from  the  lower  plateau  to  the  edge,  and  a  positive  pulse  at  the 
transition  from  the  edge  to  the  upper  plateau  (see  also  Horn  1974,  Marr  and  Hildreth  1980).  They  suggested 
that  the  image  intensities  at  the  locations  of  the  positive  and  negative  pulses  could  be  used  to  set  thresholds  to 
use  in  segmenting  the  image  into  regions. 

Several  authors  have  noted  the  relative  insensitivity  of  human  perception  to  small  intensity  gradients 
(Hcrskovits  and  liinford  1970,  Marr  1976,  Marr  and  Hildreth  1980,  McCann  ct  al.  1974).  They  have  noted 
that  the  effect  can  be  explained  by  assuming  that  the  vision  system  uses  operators  approximating  second 
derivatives.  Ihis  so-called  lateral  inhibition  effect  seems  to  be  performed  by  center  surround  operators  in 
the  retina  (sec  for  example  Richter  and  Ullman  1980).  Ihe  Laplacian  is  a  rotationally  symmetric  second 


t  A  proof  of  this  is  given  in  Scclion  3  below. 
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differential  operator,  and  an  attractive  candidate  to  perform  lateral  inhibition. 

The  use  of  the  Laplacian  for  edge  detection  was  proposed  by  Hom(1974)  in  a  study  of  the  determination 
of  lightness.  Following  Land  and  McCann(1971),  Horn  restricted  attention  to  images  of  planes  colored  with 
patches  of  uniform  reflectance  or  color.  Within  a  patch,  grey  level  variations  are  due  to  small  variations 
in  illumination,  and  they  arc  smooth  compared  to  the  abrupt  changes  between  patches.  The  conventional 
approach  to  detecting  significant  changes  in  intensity  had  been  to  note  that  the  gradient  of  the  image  is 
small  within  a  region,  but  is  infinite  across  a  reflectance  boundary  between  regions.  For  a  particular  image 
tcsselation  and  quantization  of  grey  levels,  the  gradient  is  always  finite.  It  is  usually  much  larger,  however,  at 
a  reflectance  boundary  than  it  is  within  a  region.  Hom(1974)  rejected  using  the  gradient  since  "the  first  partial 
derivatives  are  directional  and  thus  unsuitable  since  they  will  for  example  completely  eliminate  evidence  of 
edges  running  in  a  direction  parallel  to  their  direction  of  differentiation."  The  Laplacian  is  the  lowest  order 
linear  combination  of  derivatives  that  is  rotationally  symmetric.  A  reflectance  boundary  can  be  delected  by  the 
paired  positive  and  negative  peaks  on  cither  side  of  the  boundary,  and  localized  by  noting  the  position  where 
the  laplacian  crosses  zero  between  the  peaks*. 

Marr  and  Hildrcth(1980)  have  proposed  that  edges  arc  detected  in  the  human  visual  system  by  ar 
operator  that  approximates  AG,  where  A  is  the  Laplacian,  and  G  is  a  rotationally  symmetric  Gaussian.  We 
shall  show  in  the  next  section  that  the  application  of  a  rotationally  symmetric  operator,  such  as  the  Laplacian, 
to  a  rotationally  symmetric  function,  such  as  the  Gaussian,  is  itself  rotationally  symmetric.  It  follows  that  the 
Marr-Mildrcth  operator  is  rotationally  symmetric.  Morr  and  Hildreth  note  that  intensity  changes  occur  at  a 
number  of  scales  and  are  often  superimposed.  They  suggest  that  an  image  should  be  smoothed  by  a  number 
of  bandpass  filters  to  isolate  the  changes  at  a  particular  range  of  scales.  The  Gaussian  is  chosen  as  the  filter  to 
optimize  localization  of  changes  in  both  the  spatial  and  frequency  domains. 

We  noted  above  that  the  Gaussian  and  the  Laplacian  have  figured  promincmly  in  early  visual  processing 
The  Gaussian  has  mostly  been  used  to  approximate  the  point  spread  function  corresponding  to  the  blurring  of 

’  S.-c  ItmloiJil'IXi;  for  more  on  (he  dilwetion  between  (Icicciion  ,md  localisation  of  an  intensity  change. 


a  point  source.  Marr  and  Hildreth  deliberately  introduce  Gaussian  blurring.  They  further  note  that  A G  can  be 
approximated  by  a  difference  of  Gaussians,  G\  —  G^.  Nishihara  and  Larson(1981)  note  that  the  difference  of 
Gaussians  is  to  be  preferred  on  grounds  of  efficiency.  Maclcod(1972)  proposes  an  edge  detection  operator  that 
is  the  difference  of  two  Gaussians.  However,  no  analysis  of  its  performance  is  given,  and  no  indication  is  given 
that  the  operator  approximated  a  low-pass  filtered  second  derivative. 

Regarding  the  use  of  the  Laplacian,  Marr  and  Hildreth  do  not  seem  to  make  isotropy  an  explicit  con¬ 
straint  on  edge  detection.  Instead,  Hildrcth(  1980, page  13)  notes  that  "a  number  of  practical  considerations, 
which  will  be  illuminated  in  the  discussion  of  the  implementation,  suggested  that  the  . . .  operators  not  be 
directional".  Suppose  instead  that  directional  operators  are  used.  The  simplest  algorithm  for  edge  detection 
has  two  stages.  First,  the  image  is  convolved  with  the  directional  operators  in  "sufficiently  many"  directions. 
Second,  the  outputs  are  combined  to  determine  the  orientation  and  extent  of  intensity  changes.  Regarding 
{  the  first  stage,  both  Marr  and  Hildreth(1980,  page  193)  and  Hildrcth(1980,  page  40)  claim  that  the  cost  of 

convolving  the  image  with  a  "sufficient"  number  of  operators  is  excessive.  They  show  that  a  single  rotationally 
symmetric  operator  (the  I  .apiarian)  gives  precisely  the  same  results  if  a  condition  called  "linear  variation" 
holds.  Regarding  the  second  stage,  Hildrcth(1980,  page  36)  observes  that  edges  in  a  direction  close  to  that  of 
the  mask  are  elongated  in  the  direction  of  the  mask,  She  also  notes  that  operators  at  several  orientations  give 
significant  responses  to  any  given  edge,  and  that  combining  the  responses  is  non-trivial. 

There  are  two  essentially  different  issues  here  that  need  to  be  clearly  separated.  Intensity  changes  first 
have  to  be  detected  and  then  localized  as  a  set  of  "feature  points"  marking  the  position  of  the  change  in  the 
image,  and  characteristics  of  the  corresponding  edge.  detection  of  feature  points  is  inherently  isotropic, 
as  Hom(1974)  noted.  The  feature  points  have  then  to  be  combined  to  produce  descriptions  of  ed/e  segments. 
Edge  segments  arc  clearly  directional,  indeed  a  central  problem  concerns  the  determination  of  the  direction  of 
an  edge  in  an  image.  The  compulation  of  rich  descriptions  of  edge  segments  is,  as  Hildreth  notes,  not  at  all 
easy.  Marr’s(1976)  original  Primal  Sketch  work  was  almost  entirely  concerned  with  it.  Uinidrdf  1981)  discusses 

1 

the  application  of  directional  operators  to  compute  the  directionality  of  an  edge. 


The  Gaussian  and  Laplacian  arc  not  the  only  rotationally  symmetric  operators  that  have  been  proposed 
in  computer  vision.  Prcwitt(1970,  p.  107)  observes  that  "derivatives  of  all  orders  can  be  used  to  form  isotropic 
nonlinear  differential  operators,  provided  that  derivatives  of  odd  order  appear  only  in  even  functions.  The 
simplest  of  these ...  is  the  squared  gradient",  namely  V  •  V.  where  V  is  the  column  vector 

■&' 

A 

Earlier  in  the  same  article,  Prewitt(1970,  p.  85)  suggests  that  "the  Hankcl  transformation  enters  naturally 
in  the  analysis  of  systems  with  isotropic  point  spread  functions  and  greatly  facilitates  restoration."  The  sugges¬ 
tion  does  not  appear  to  have  been  investigated  in  computer  vision. 

We  noted  earlier  that  an  important  aspect  of  modelling  perception  is  the  isolation  of  constraints  which 
capture  some  facet  of  smoothness.  Horn  and  Schunck(1981)  consider  the  determination  of  optical  flow  fields 
and  note  that  "if  every  point  of  the  brightness  pattern  can  move  independently,  there  is  little  hope  of  recover 
ing  the  velocities".  One  way  to  express  the  additional  constraint  of  smoothness  is  to  minimize  the  integral  of 
Jie  performance  index 


S(u,  v )  =  (u*  -f  «})  4  («2  +  «J). 

where  u  and  u  are  the  x  and  y  components  of  the  optical  flow,  and  subscripts  denote  partial  differentiation. 
We  shall  show  in  the  next  section  that  this  operator  is  rotationally  symmetric.  In  many  simple  situations  the 
smoothness  constraint  is  significantly  wrong  only  at  occluding  boundaries. 

Wc  conclude  this  review  of  the  use  of  rotationally  symmetric  operators  in  vision  with  Grimson’s(1981) 
woik  on  surface  interpolation.  As  t'<  will  be  the  focus  of  Section  5.  our  remarks  will  be  brief.  The  Marr-Poggio 
theory  ol  human  stereo  vision  yields  the  di*  parity  (scaled  depth)  at  matched  edge  points  that  arc  computed 
by  the  Marr-Hildrcth  approach  described  above.  The  disparity  map  is  as  sparse  as  the  set  of  matched  edge 
points,  whereas  human  perception  is  of  smooth  surfaces  passing  through  the  given  disparity  points.  Grimson 
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(1981)  interpolates  a  smooth  surface  from  the  given  set  of  edge  points  by  a  local  parallel  algorithm  that  applies 
a  rotationally  symmetric  operator  to  minimize  die  quadratic  variation  introduced  above. 

3.  Conditions  for  rotational  symmetry 

A  function  *-*  R  is  rotationally  symmetric  if  its  polar  form  is  only  dependent  on  radial  distance 
r  =  (i2  -f-  y2)4  and  not  on  direction  <j>  =  tan-1  ~ .  Clearly,  a  function  is  rotationally  symmetric  if  and  only 
if  it  can  be  represented  as  a  function  of  (*2  -}-  y2)4,  An  alternative  definition  can  be  given  that  is  often  more 
convenient  for  functions,  and  that  can  be  generalized  to  operators.  A  ftmction  is  rotationally  symmetric  if  and 
only  if  it  yields  the  same  value  under  an  arbitrary  rotation  of  coordinates. 

An  anticlockwise  rotation  from  one  set  of  image  coordinates  (*,y)  to  another  (X,  T)  is  effected  by  a 
rotation  matrix: 


ain^' 
cos  <f>. 


(0) 


For  convenience,  we  shall  denote  cos^  by  c  and  sin  0  by  «.  To  simplify  notation,  we  shall  not  make 
explicit  the  dependence  of  the  rotation  matrix  R  on  the  angle  <t>.  A  function  /  is  rotationally  symmetric  if  and 
only  if  the  untransformed  version  f(x,  y)  is  equal  ^  the  transformed  version  f(X,  Y).  We  shall  occasionally 
find  it  useful  to  borrow  the  mathematical  shorthand  that  equates  a  function  f(X,  Y)  with  a  function  of  a  single 
vector  argument  f{R\x,  y]r). 

Example  /.  The  function  f\{x,  y)  =  (*2  -f  V2)  is  rotational^  symmetric: 


MX,  Y)  —  ((arc  +  y«)2  +  (yc  —  m)3) 

=  (*2  +  vJ) 

=  f\[x,  y). 


Example  2.  The  function  /.(z,  y)  =  xy  is  not  rotationally  symmetric: 


8 


MX,  Y)  *  (*c  +  ya)(yc  —  xt) 

—  ** 

*  xy  coi  2<fi  +  — j —  ain  2^, 

and  so  h(X>  Y)  —  Mx>  v)  only  when  ^  =  0  or  ^  =*  ir. 

We  can  extend  the  definition  of  rotational  symmetry  to  operators 

0:(B2  w  *)  *-♦  (Ra  ►-  R). 

An  operator  0  is  rotationally  symmetric  if  0(/)  is  a  rotationally  symmetric  function,  for  all  functions 
/:£*  h*  *. 

Example  3.  The  function  produced  by  the  operator  Oj,  defined  by 

Oi(/)(x,y)»e^*-v) 

is  rotationally  symmetric  if  and  only  if  /  is.  In  general  then,  the  operator  Oj  is  not  rotationally  symmetric. 
However,  the  Gaussian  is  a  rotationally  symmetric  operator  as  it  combines  examples  1  and  3. 

Most  of  the  operators  of  interest  in  computer  vision  are  combinations  of  the  first  and  second  directional 
derivatives  and  We  need  to  determine  the  effect  of  a  coordinate  rotation  on  these 

directional  derivatives.  By  the  chain  rule, 

d_  =dXd_  dYd_ 

dx  dx  dX  dx  dY 

d  d 

CdX  'dY 

Similarly, 

d  d  ,  d 

dy  ex  ^  dY 


•«fcvT' 
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[X  f*]V{X,Y)  =  [x  m]V(,,„). 


By  equation  (2), 


(X  =  [X 

and  so  the  linear  differential  operator  is  rotationally  symmetric  if  and  only  if 


(X  M]  =  [X  h]R, 

so  that  [,  /x]  :.s  an  eigenvector  of  R.  The  eigenvalues  of  R  are  c  -f  is  and  c  —  is.  So  there  are  no  real 

eigenvectoi '  unless  <f>  is  a  multiple  of  ir.  Since  the  condition  is  not  satisfied  for  all  <£,  no  linear  combination  is 
rotationally  symmetric.| 

The  same  style  of  analysis  can  be  applied  to  other  combinations  of  first  derivatives  such  as  the  operator 


02(/)  = 


U~fv 

/*  +  /» 


|  It  is  easy  to  show  that  02(x,r)  is  not  equal  to  02(Xiy),  for  example  when  ^  =  f . 

! 

In  section  2.  we  referred  to  an  operator  proposed  by  Prewitt(1970),  namely 


that  is.  the  vector  dot  product 


More  generally,  wc  often  consider  quadratic  differential  expressions  such  as 


. a*cas‘-- '  a 


rx 

U'  SJ 

Such  an  expression  is  called  a  quadratic  form  if  the  matrix  is  symmetric,  that  is  n-v.  By  equation  1, 


so  that 


V(x,y)  = 


if  and  only  if 


where  R  is  an  arbitrary  rotation  matrix,  and 


RtMR  =  M, 


rX  Ml 

M  = 

\y  eJ 

Since  the  transpose  RT  of  a  rotation  matrix  R  is  the  inverse  of  R,  a  quadratic  form  is  rotationally 
symmetric  if  and  only  if  the  corresponding  matrix  M  commutes  with  all  rotation  matrices.  We  will  refer  to 
matrices  M  having  this  property  as  being  rotationally  symmetric. 

licmma  I.  A  2  by  2  matrix  is  rotationally  symmetric  if  and  only  if  it  has  the  form 


M  = 


X  M 
—H  X. 


Proof.  We  require  RM  =  MR  for  all  rotation  matrices  R,  that  is 


.Mil 
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'  ) 

c  — «irx  /ii  rx  /iirc  — j 

*  c  JU  d  U  dl*  c . 

Expanding,  and  equating  terms,  this  holds  if  and  only  if 

x  =  e- 

Alternatively,  only  the  operations  of  scaling  by  a  constant  k  and  multiplication  by  a  rotation  matrix  Bl 
commute  with  all  rotation  matrices  in  two  dimensions  So  M  —  kR'  for  some  scale  factor  k  and  some  rotation 
matrix  R'i 

Proposition  2.  Up  to  scaling,  the  only  rotationally  symmetric  quadratic  form  in  £  and  ^  is  V(I>v)  •  V(I>W). 

Proof.  A  quadratic  form  in  and  ^  has  the  form 


To  be  rotationally  symmetric,  as  well  as  symmetric  (so  that  it  is  a  quadratic  form).  Lemma  1  implies  that 

x  =  e 

^  =  0. 

It  follows  that  the  matrix  in  equation  (3)  is  X/jj 

The  operator  f\- 1-  /J  is  commonly  used  as  a  measure  of  the  contrast  across  an  intensity  change.  Notice 
that  other  popular  measures  of  the  contrast,  such  as  {/*  +  /„)2,  (fz  —  /y)2,  or  ||/?||  -f  ||/J|  arc  not  rotationally 
symmetric,  and  therefore  respond  differently  to  edges  in  different  directions  (see  Rosenfcld  and  Kak  1976, 
p279) 

Wc  now  consider  linear  and  quadratic  forms  in  and  ^ .  It  is  convenient  to  not  assume 

=  f°r  the  developments  that  follow. 

The  first  task  is  to  find  a  matrix  /?*  so  that 


O 
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(4) 


The  (*,  j)  element  of  the  matrix  R  t  will  be  denoted  by  r,  j.  Applying  the  chain  rule  as  before,  but  this 
time  to  relate  the  second  derivatives  in  (X,  Y)  to  those  in  (*,  y),  we  find  that  the  four  by  four  matrix  R*  can  be 
written  in  the  form 


t 


r* 


rnilr  rnRT 

f|2-Rr  f22 RT 


Definition  1.  (ben  Israel  and  Greville  1974,  page  41)LetA  —  [o,;]  and  B  —  [6tJ]  be  m  by  m  and  n  by  n 
matrices  respectively.  The  mn  by  mn  matrix  A  ®  B,  called  the  Kronecker  product  of  A  and  B,  is  defined  by 
multiplying  each  element  a(t,  j)  of  A  by  the  matrix  B,  to  form  die  block  matrix 


on  B 

oi»B  . . . 

«1  ttJB 

021 B 

022$ 

02mP 

flm\B 

OmlB  .  . . 

OmrnB 

With  this  notation, 

R'  =  Rt®Rt, 


(«) 


so  that 


a 


1 


'  Recall  the  definition  of  the  matrix  R  from  equation  (0) 
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Note  that  the  elements  of  /4  <g>  J9  are  naturally  indexed  by  4*tuples: 


(7) 


{A<g)B}ijkl  =  Oijiu. 

Wc  state  without  proof  a  number  of  simple  properties  of  the  ®  operation.  They  are  essentially 
straightforward  consequenc  es  of  the  properties  of  ordinary  multiplication,  and  are  stated  without  proof. 
Proposition  3 

(t)  [A  ®2?)t  =  At  ®Bt 
(«)  [A  ®  B)~*  =  A~l  <S)B~l 
Jiii)  (A®B)<g)C  =  A<8>(B<g)C) 

For  the  remainder  of  the  paper,  we  restrict  attention  to  the  application  of  ®  to  R  and  its  transpose. 
Proposition  4.  The  rotationally  symmetric  linear  combinations  of  sSy*  fife-  and  $  are  linear 
combinations  of  the  Laplacian  A  =  jfp  +  and  the  smoothness  measure  gj§j  —  g^j. 

Proof.  Let  the  linear  combination  be 

[*  M  v  l\ 

Following  the  proof  of  Proposition  1,  the  condition  for  rotational  symmetry  is 


I*  M  v  t\RT(g>Rr=[\  (i  u  $], 


IS 


for  all  rotation  matrices  R  and  the  corresponding  rotation  angle  Expanding  RT  ®  RT  by  equation  (7),  we 
find 


[X  M  €] 


c2  —bc  —Be  a2 


BC 


BC 


c2  -a2 


-i* 


BC 


C‘  —BC 


—  I*  M  V  6]» 


so  that 


[X  /i  i/ 


It  follows  that 


L2 

BC 

BC  C2  J 

—BC 

— »c  a2 

BC 

-a2 

— a2  — bc 

BC 

-a2 

1 

1 

S 

.a2 

BC 

bc  — a2J 

r- 

2 a2  - 

2ac  — 2ac 

=  [0  0  0  0]. 


[X  —  C  M  +  *'  0  0) 


2ac  —2a2 
0  0 

L  0  0 

The  determinant  of  the  upper  left  2  by  2  submatrix  is 


-2a2 

0 

0 


-2ac 

0 

0  . 


=  [0  0  0  0). 


(4a*  -f  4a2c2)  =  4aa. 

Since  this  is  not  zero  for  all  angles  X  —  £  and  n  v  are  both  zero.  A  basis  for  the  infinite  set  of 
linear  combinations  satisfying  these  conditions  is  provided  by  setting  X  and  n  equal  to  one,  which  proves  the 
Proposition.  | 

Before  turning  to  quadratic  forms,  analogous  to  Proposition  (2),  we  define  a  projection  operator  on 
Rr  (gi  R1  that  makes  explicit  the  assumption  —  fv*- 
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Definition  2.  Let  D  =  \di}\  be  a  4  by  4  matrix.  The  projection  of  D  is  the  3  by  3  matrix  D* : 

<* ii  (<*12  4  <*u)  <*14 

(<*21  4-  <*3l)  (<*22  -f  <*32  -4  <*23  4  <*33)  (<*24  +  <*3*) 

<*41  (<*42  +  ^43)  <*44 

That  is,  the  second  and  third  columns  as  well  as  the  second  and  third  rows  are  combined  by  addition. 
Proposition  5. 

[abb  c]2?[o  b  b  c]T 


is  equivalent  to 


(o  6  cp*(o  6  c]T, 

where  D*  is  the  projection  of  D. 

TTic  proof  is  by  equating  terms,  and  is  omitted.  Wc  now  give  the  main  result  of  this  section. 

Proposition  6.  The  rotationally  symmetric  quadratic  forms  in  g^,  g^gj,  and  form  a  vector 
space.  If  —  gjjgj,  the  matrices  associated  with  the  rotationally  symmetric  quadratic  forms  project  to  3  by 
3  matrices  of  the  form 


Q-f  0  0  0 

0  2a  0 

0  0  a40 

It  follows  that  the  rotationally  symmetric  quadratic  forms  that  satisfy  form  a  vector  space 

that  has  the  quadratic  variation,  (^)2  4  2(g^)2  4  (^)2,  and  the  square  Uplacian,  (^  4  J^s)2.  as  a 
basis. 

Proof.  Since  the  matrix  in  a  quadratic  form  is  defined  to  be  symmetric,  a  quadratic  form  in  g^, 
0%  ■  and  can  be  written 


• .. .  s.  "  rw.  -  - "  r-i. 


'.4? 


$  Sy&  $\ 

J  ^6  C \ 


C\  flyM 


where  A  and  C  are  symmetric  2  by  2  matrices,  and  B  is  2  by  2.  As  usual,  the  quauratic  form  is  rotationally 


symmetric  if  and  only  if 


rt ® rt\ A  B]  =  rA 

L BT  Cl  BT  C\ 


where  R  is  an  arbitrary  rotation  matrix.  It  follows  that 


cRt  aR 


-aRr  cRtj 


nr  4  b]  [a  sir  a 

tabt  c\  Ur  cl.-, 


cRT  aR7 


aRT  cRTl 


and  hence  that 


'  cRtA  4-  sRrBT  cRtB  4-  aRTC  cART  -  aBRT  aARr  +  cBRT 
—aRrA  4-  cRtBt  —bRtB  4-  cRtC  cBtRt  -  aCRT  aBTRT  4  cCRT 

Equating  submatriccs,  we  find  that  for  all  rotation  angles  4> 


c(RT A  -  ARt)  4  a(RTBT  +  BRt)  =  0, 
c(RTC  -  CRr)  -  a(BTRT  4  RTB)  =  0, 
a(CRT  -  RtA)  4  c(RTBT  -  BtRt)  =  0, 
a{RrC  -  AR1')  4  c(RtB  -  BRT)  =  0. 

Consider  equation  (10)  or  (11)  when  <j>  —  f .  Equating  terms,  we  find  that 

on  =  C22 

022  =  cji 

0|2  =  —  C2l 
021  a  — C)2. 


(12) 


Similarly,  equation  (8)  or  (9)  when  ^  **  \  yields 


&11  4"  &22  =»  0. 

Expanding  equation  (8)  for  general  <p  yields 


(13, 


bu  4-ou  —  0, 

622  —  aji  =  0, 

&21  +  bi2  +  022  —  flu  —  0. 

Combining  equations  (12)  through  (16)  we  find  that  in  order  to  be  rotationally  symmetric,  the  matrix 


(14] 

(15] 
(1«] 


has  the  form 


A  matrix  of  this  form  projects  to 


'<*  +  0  7 

7  a  —  S 
7  6 

.  P  7 


—1  h' 

6  7 

a  — S  — 7 

Q  +  P. 


o  +  /J  0  0 

0  2a  0 
.  0  0  a  +  0 

where  a  =  bi2  —  on  and  0  =  bl2.  It  is  easy  to  show  that  linear  combinations  of  matrices  of  this  form  arc 
ol  the  same  form,  so  that  the  rotationally  symmetric  quadratic  forms  constitute  a  vector  space.  Clearly,  the 
square  I  .apiarian  and  the  quadratic  variation,  corresponding  to  the  cases  c  =  1,0  =  0  and  a  —  0,0  =  1 
respectively,  form  a  basis* 
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We  show  that  the  measure  of  smo  Jiness  of  optical  flow  proposed  by  Horn  and  Schunck(1981)  is  rota* 
tionally  symmetric.  Recall  from  section  2  that  the  measure  is  defined  by  the  operator 


s(«,u)  =  (U» +„;)+(»;+»}). 

We  extend  the  Kronccker  product  operator  ®  to  vectors,  and  then  show  how  to  define  S(u,  v)  in  terms 
of  vector  Kronecker  products. 

Definition  3.  (a)  Let  a  =  (oi. .  .a*,]  and  fe  =  [bt .  .bn]  be  vectors.  The  Kronecker  product  of  a  and  fe  is  the 
mn  dimensional  vector  [aj6i. .  aibn  a2fci, .  .Ombn], 

(b)  By  extension,  ifQ  =  (0). .  .0mj  is  a  vector  of  operators  and/  =  |/i.../n]  is  a  vector  of  functions,  the 
Kronecker  product  of  <1  and  £  is  the  mn  dimensional  vector  of  functions 


[01(/1)...01{/n)...0m(/n)]. 

The  components  u  and  v  of  optical  flow  are  functions  of  *,  y,  and  t.  Recall  th?t  V(*iy)  =  \§^  ^]r. 

According  to  definition  \ 


so  that 


du  dv 

(9y  dx 


S(u,  v)  =  (v{,iW)  ®  [u  v]T)  •  (V  ,,y)  ®  [u  v]T ). 

If  the  coordmate  frame  is  rotated  through  4>  by  the  mafix  R,  the  optical  flow  components  become  R\u  v]T. 
The  Horn-Schunck  measure  is  ^nationally  symmetric  if  and  only  if 


{R®R)t{R®R)=U, 
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where  U  Is  tiie  4  by  4  identity  matrix.  The  rotational  symmetry  is  a  simple  consequence  of  Proposition  3. 
A  rotationally  symmetric  operator  has  the  general  form 


and  its  application  to  a  rotationally  symmetric  function /(x,  y)  has  the  form 

Ci»’ 


0(*,v)(/(*»  v))‘ 

To  see  that  this  is  rotationally  symmetric,  we  rotate  the  coordinate  frame  to  (X,  T)  by  a  matrix  R  as  before. 
Since  0  and  /  arc  rotationally  symmetric,  all  the  occurences  of  R  (including  its  Kronccker  square,  cube,  and  so 
on)  introduced  by  the  frame  change  can  be  deleted.  It  follows  that  the  application  of  a  rotationally  symmetric 
operator  to  a  rotationally  symmetric  function  is  itself  rotationally  symmetric.  In  particular,  the  A(C7)  filters  of 
the  Marr-Hildrcth  tneory  of  edge  detection  are  rotationally  symmetric. 

4.  Vision  as  a  conservative  process 

The  second  theme  of  this  paper  is  that  a  number  of  vision  modules  construct  the  most  conservative  inter¬ 
pretation  that  is  consistent  with  the  given  data,  and  that  is  subject  to  an  appropriate  set  of  suitably  formulated 
constraints.  A  major  concern  of  Computer  Vision  has  always  been  the  isolation  of  constraints  that  enable  the 
interpretation  of  an  image.  Constraints  embody  observations  about  the  way  the  world  is,  at  least,  most  of  the 
time.  Although  such  observations  can  be  as  specific  as  cataloging  familiar  figures  and  shapes,  it  has  proved 
more  fruitful  to  first  uncover  constraints  that  correspond  to  general  observations  that  are  widely  applicable. 
C  onst  taints  arc  used  together  with  the  data  computed  from  the  image  to  construct  an  interpretation.  The 
representations  of  the  information  from  the  image  and  tire  constraints  determine,  and  arc  determined  by, 
the  interpretation  process.  For  example,  early  blocks  world  programs  represented  constraints  as  catalogs  of 
labellings,  an  apptoach  that  led  naturally  to  search  processes  for  interpretation  (Clowes  1971,  Kanadc  1981). 
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As  Computer  Vision  has  turned  its  attention  to  images  of  the  natural  world,  constraints  have  concerned 
the  smoothness  of  surfaces  and  movement  The  relationship  to  boundary  value  problems  of  physics  and 
mathematics  suggests  itself.  The  information  computed  from  the  image  sets  the  boundary  conditions,  and  the 
constraints  determine  which  (and  whether  a)  solution  to  the  boundary  value  problem  is  found.  Hom(1974) 
solved  an  instance  of  Poisson's  problem  using  Green's  functions  to  determine  the  lightness  of  an  image. 

Following  a  different  approach,  Ullman(1979a)  studied  the  perception  of  apparent  motion  generated 
by  two  successive  frames  consisting  of  isolated  dots  of  equal  intensity  moving  independently  of  each  other. 
Without  constraint,  all  possible  pairings,  or  "correspondences",  of  dots  in  the  first  frame  with  dots  in  the 
second  are  equally  likeiy.  Ullman  defined  the  "most  likely"  correspondence  to  be  the  one  that  minimized  the 
sum 


t<«;Sn 

where  n  is  the  number  of  dots  in  the  first  frame,  m  is  the  number  of  dots  in  the  second  frame,  and  is  one  if 
the  ith  dot  of  the  first  frame  P,  is  paired  with  the  ;th  dot  of  the  second  frame  Qj,  else  zero.  The  weight  fty  is 
the  cost"  of  pairing  P,  with  Qj,  and  might,  for  example,  be  related  to  the  image  distance  between  the  paired 
points.  The  problem  of  finding  the  minimal  correspondence  is  considered  in  terms  of  integer  programming.  If 
correspondences  are  assumed  to  be  covering  mappings,  the  following  linear  constraints  apply  to  the  *y: 

Vi,l<i<n  £ 


and 


V?*  1  <  J  <  m  £  Xij  >  l. 

I<i<n 


Ullman  restricted  the  set  of  Qj  that  can  be  paired  with  to  those  whose  positions  were  close  to  Following 
Arrow,  Hurwicz,  and  Uzawa(19S8),  he  set  up  the  iterative  scheme 


“!+l  -  E  *!i-> 

*5+'  =  £  ~  1 

The  approach  can  be  extended  to  mapi  mgs  that  are  not  one-one,  as  well  as  to  continous  motion.  A  major 
problem  with  the  approach  is  guaranteeing  the  convergence  of  the  algorithm.  This  is  determined  largely  by  the 
properties  of  the  costs  g,j,  but  this  was  not  investigated,  aside  from  a  comment  on  the  empirical  determination 
of  the  qtJ  (see  also  Ullman  1979b). 

One  limitation  of  Ullman ’s  approach  is  that  it  is  restricted  to  minimizing  a  known  linear  objective  func¬ 
tion  that  is  subject  to  linear  constraints.  The  method  can  be  extended  to  constrained  nonlinear  programming 
in  which  the  goal  is  to  minimize  a  known  function  /(f)  subject  to  a  set  of  equality  and  inequality  constraints 
of  the  form  g,(f)  <  0.  In  general,  however,  criteria  based  on  other  than  intuition  need  to  be  found  for 
selecting  the  function  /  to  be  minimized.  To  do  this,  one  can  apply  the  calculus  of  variations  (see  for  example 
Courant  and  Hilbert  1953,  chapter  IV).  The  familiar  differential  shows  how  to  find  a  real  valued  parameter 
that  minimizes  some  function.  The  calculus  of  variation  extends  the  differential  calculus  by  showing  how  one 
can  determine  a  Junction/*,  which  is  subject  to  a  given  set  of  boundary  conditions,  and  minimizes  the  integral 

*^(/)  ==  f  JG^X’  V'  fxx,  fxy,  fvv)dxdy  (17) 

over  a  given  region  of  integration  G\  The  function  F  is  called  a  "performance  index"  and  generalizes  the 
notion  of  cost  function  associated  with  linear  and  nonlinear  programming.  In  the  next  section  we  shall  con¬ 
sider  the  choice  of  a  performance  index  for  interpolating  smooth  surfaces  from  one-ditncnsional  boundary 
conditions. _ 

1  l  or  simplicity  of  presentation,  we  rcstiict  attention  to  functions  /  of  one  or  two  variables  x ,y. 
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Associated  with  a  variational  problem  of  the  form  (17)  is  the  Euler  equation,  which  provides  a  necessary, 
though  by  no  means  sufficient,  condition  which  a  function  /  must  satisfy  if  it  is  to  minimize  the  variational 
integral  ?(/).  For  the  particular  variational  problem  given  in  equation  (17),  the  Euler  equation  is 

Ff  ~  dxFf‘  ~  dyF/y  dx *F/mm  dxdyFf*v  "**  d\?Ffvy  “ 

In  the  case  that  there  is  only  a  single  dependent  variable  x,  the  partial  derivatives  are  total  and  the  Euler 
equation  becomes 

+  (1») 

There  are  two  important  considerations  associated  with  the  use  of  the  calculus  of  variations.  First,  unlike 
the  differential  calculus,  the  existence  of  an  extremum  /*  of  the  integral  given  in  equation  (17)  cannot  be  taken 
for  granted.  Courant  and  Hilbert(1953,  p.  173)  note  that  "a  characteristic  difficulty  of  the  calculus  of  variations 
is  that  problems  which  can  be  meaningfully  formulated  may  not  have  solutions".  Conditions  for  the  existence 
of  a  minimum  have  recently  been  proposed  by  Grimson(1981)  and  will  be  discussed  in  the  next  section. 

Second,  associated  with  any  variational  problem  is  a  set  of  natural  boundary  conditions  which  imposes  a 
necessary  condition  on  any  feasible  solution  to  the  Euler  equation  at  the  boundary.  Courant  and  Hilbert(1953, 
p.  211)  note  that  "in  general,  we  can,  by  adding  boundary  terms  or  boundary  integrals  essentially  modify 
the  natural  boundary  conditions  without  altering  the  Euler  equations".  Determining  the  "most  conservative" 
solution  means  finding  a  performance  index  that  guarantees  the  existence  of  an  extremum  function  /*  and 
provides  the  tightest  set  of  natural  boundary  conditions  that  are  consistent  with  the  given  data. 

The  calculus  of  variations  has  recently  been  applied  by  a  number  of  authors  to  interpolate  plane  and 
space  curves  and  surfaces.  We  review  the  applications  in  that  order.  First,  Horn(1981)  has  recently  determined 
the  curve  which  passes  through  two  specified  points  with  specified  orientation  while  minimizing 


where  n  is  the  curvature  and  *  is  the  arc  length.  This  is  the  true  shape  of  a  spline  used  in  "lofting"  (Faux  and 
Pratt  1979, p.  228).  In  a  thin  beam,  curvature  is  proportional  to  the  bending  moment.  The  total  elastic  energy 
stored  in  a  thin  beam  is  therefore  proportional  to  the  integral  of  the  square  of  the  curvature.  Sipce  the  shape 
taken  on  by  a  thin  beam  is  the  one  which  minimizes  the  internal  strain  energy,  the  curve  that  solves  equation 
(20)  is  called  the  "curve  of  least  energy".  The  variational  problem  is  to  minimize 

f  --**  idi- 
1  u+/a! 

This  has  the  form  of  equation  (17).  Horn(1981,  page  19)  shows  that  the  Euler  equation  is 

— ck  —  y  cos  V>i 

where  t/>  is  the  angle  between  the  tangent  to  the  curve  and  the  axis  of  symmetry.  The  solution  to  this 
differential  equation  is  an  incomplete  elliptic  integral  of  the  first  kind.  Brady,  Grimson,  and  Langridge(1980) 
consider  a  "small  angle"  approximation  to  the  curve  of  least  energy,  in  which  first  derivatives  can  be  ignored. 
The  performance  index  that  they  use  is  f\x,  for  reasons  that  will  become  evident  in  the  next  section.  They  find 
that  in  that  case  the  solution  is  a  cubic.  Hom(1981,page  2)  notes  that  the  fact  that  a  curve  has  near  minimum 
energy  does  not  mean  that  it  lies  close  to  the  curve  of  minimum  energy.  Note  that  the  existence  of  the  curve 
of  least  energy  is  guaranteed  as  Horn  has  derived  an  analytical  formula  for  it.  Approximations  to  it,  such  as 
Brady,  Grimson,  and  Langridgc’s  arc  similarly  guaranteed  to  exist. 

Bai  iow  and  Tencnbaum(1981)  investigate  the  problem  of  interpreting  a  line  as  the  image  of  a  space  curve 
dial  is  an  occluding  boundary.  'ITiey  observe  that  the  problem  has  two  parts:  (i)  determining  the  tangent 
vector  ( at  each  point  on  the  space  curve,  and  (ii)  determining  the  surface  normal  at  each  point,  given  that  it  is 
constrained  to  be  orthogonal  to  the  tangent. 


v...  .Vr-W  ,  IWWIlliHH  rt  WS^jjl^gfegfe81*  - 
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They  suggest  minimizing  a  performance  index  F  that  is  a  function  of  the  curvature  k  and  the  torsion  r 
(possibly  together  with  their  derivatives),  and  expresses  a  suitable  notion  of  "smoothness".  They  first  consider 
uniformity  of  curvature  as  a  measure  of  smoothness,  that  is  F  =  ^  =  /cj,  where  a  measures  distance  along 
the  space  curve.  They  reject  this  measure  on  the  grounds  that  n\  can  be  made  arbitrarily  small  by  "stretching 
out  the  space  curve  so  that  it  approaches  a  twisting  straight  line".  To  overcome  this  difficulty,  they  propose 
that  the  space  curve  should  also  be  "as  planar  as  possible  or,  more  precisely,  that  the  integral  of  its  torsion 
should  be  minimized". 

Barrow  and  Tenenbaum  finally  suggest  finding  the  space  curve  that  projects  to  the  given  image  line  and 
minimizes  the  performance  index  [^p]?,  where  fc  is  the  binormal.  They  report  that  an  algorithm  based  on 
their  analysis  produced  the  "correct  3-D  interpretations  for  simple  and  closed  curves,  such  as  an  ellipse,  which 
was  interpreted  as  a  circle".  However,  they  note  that  the  rate  of  convergence  was  slow  and  dependent  on  the 
initial  data.  No  consideration  is  given  to  the  Euler  equations,  to  the  existence  of  an  extremum  given  a  line 
drawing  (x(«),  y(fi)},  or  to  the  natural  boundary  conditions  associated  with  the  performance  index  [^p]2. 
Empirical  evidence  that  the  method  works  on  a  number  of  simple  test  cases  is  encouraging;  but  there  is  no 
analysis  of  the  scope  of  the  method. 

hi  the  same  paper,  Barrow  and  Tcnenbaum(1981)  consider  the  interpolation  of  a  smooth  surface  from 
dep-t  and  local  surface  orientation  values  at  all  points  along  the  surface  boundary.  Their  approach  is  to 
"seek  a  technique  that  yields  exact  reconstructions  for  the  special  symmetric  cases  of  spherical  and  cylindrical 
surfaces,  as  well  as  intuitively  reasonable  reconstructions  for  other  smooth  surfaces."  (Barrow  and  Tenenbaum 
1981).  They  observe  that  if  n  is  the  surface  normal  of  a  cylinder,  then  the  x  and  y  components  of  the  normal 
nr  and  n,y  are  linear  functions  of  x  and  y,  so  long  as  the  axis  of  the  cylinder  lies  in  the  x  —  y  plane.  This 
observation  forms  the  basis  of  an  algorithm  to  estimate  the  surface  normal  by  least  squares  fitting  of  the 
parameters  of  the  partial  derivatives  of  the  normal.  As  before,  no  analysis  is  given  of  the  Euler  equation,  the 
natural  boundary  conditions,  nor  the  convergence  of  their  algorithm  for  different  types  of  surface. 
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5.  A  performance  index  for  surface  interpolation. 

In  the  review  of  the  application  of  the  calculus  of  variations  to  visual  perception  in  the  previous  section 
we  drew  attention  to  three  important  considerations.  First,  .h ;  Euler  equations  provide  a  necessary  condition 
on  possible  extremal  functions.  Second,  the  existence  of  an  extremum  cannot  be  taken  for  granted,  even  when 
the  minimization  problem  seems  plausible  on  some  grounds.  Third,  the  natural  boundary  conditions  impose 
a  necessary  condition  on  any  feasible  solution  to  the  Euler  equation  at  the  boundary.  The  most  thorough 
analysis  of  the  second  of  these  problems  in  Computer  Vision,  framed  in  the  context  of  surface  interpolation,  is 
due  to  Grimson(1981),  who  proves  the  following  theorem. 

Theorem  (Grimson,  sec  Rudin(1973))  Suppose  there  exists  a  complete  semi-norm  F  on  a  space  of  func¬ 
tions  5,  and  that  F  satisfies  the  parallelogram  law.  Then,  every  non-empty  closed  convex  set  6  Cl  7  contains  a 
unique  clement/*  of  minimal  normF(/*),  up  to  possibly  an  element  of  the  null  space  of  F. 

A  semi-norm  F  is  a  function  V  h-»  Si+  from  a  vector  space  V  to  the  positive  real  numbers  that  satisfies 

F{v  +  w)<F{v)  +  F[w) 

F(qv)  =  ja|F(v). 

Informally,  a  semi-norm  is  a  generalization  of  the  Euclidean  metric,  and  provides  a  measure  of  a  vector.  The 
second  condition  generalizes  the  triangle  inequality,  for  example.  The  null  space  of  the  semi-norm  F  consists 
of  all  those  vectors  Uu  that  map  to  zero.  Since 


F{v  +  to)  =  F[v), 

any  clement  of  the  null  space  can  be  added  to  a  vector  of  minimal  norm  to  yield  another  vector  of  minimal 
norm.  Hence  the  qualifying  phrase  "unique  ...  up  to  possibly  an  element  of  the  null  space  of F".  The 
parallelogram  law  states  that 

[F(v  +  u>)]2  +  [F(v  -  ti,)]2  =  2[F(i>)]2  +  2[F(tt/)]2, 
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for  all  vectors  t>,  w.  Finally,  the  semi-norm  is  complete  if  all  Cauchy  sequences  converge.  As  is  well  known, 
the  elements  of  vector  spaces  can  be  functions.  This  enables  Grimson  to  prove  the  following  Corollary,  that 
guarantees  the  existence  of  an  extremum  function  in  calculus  of  variations  "most  conservative"  interpolation 
problems. 

Corollary  (Grimson  191(1).  Let  the  set  of  known  points  be  {(**,  y„  *,)  1 1  <  *  <  n}.  Let  9  be  a  vector 
space  of  possible  functions  on  R2  and  let  8  be  the  subset  of  9  that  interpolates  the  known  data.  That  is,  for  all 
functions  /e8,/(xj,  y.)  =  a,.  Let  F  be  a  complete  semi-norm  on  8  that  satisfies  the  parallelogram  law.  Then 
there  exists  a  unique  (up  to  the  null  space  of  F)  function  f  that  interpolates  the  data  and  has  minimal  norm. 
In  particular,  if  F  is  a  performance  index  then  there  is  a  function  /*  that  minimizes  the  integral 

9in  =  Jf 

In  short,  if  the  conditions  of  the  Corollary  arc  fulfilled,  the  existence  of  a  "most  conservative"  surface  that 
meets  the  boundary  conditions  is  guaranteed.  As  we  shall  see,  the  condition  of  being  a  semi-norm  is  the  most 
restrictive  required  of  the  performance  index.  The  conditions  are  sufficient  to  guarantee  the  existence  of  a 
minimum,  but  they  are  not  necessary.  For  example,  k 2  is  not  a  seminorm*;  nevertheless  Hom’s(1981)  analysis 
shows  that  there  is  a  unique  minimum.  It  is  far  from  clear  whether  Barrow  and  Tcnenbaum’s(  1981)  analyses  of 
curve  and  surface  interpolation  have  a  guaranteed  minimum  in  all  cases. 

Grimson  notes  that  several  intuitively  plausible  performance  indices  are  not  semi-norms.  For  example, 
the  two  most  popular  measures  of  curvature  are  not.  Suppose  that  «i  and  are  the  principal  curvatures  of 
a  surface(Faux  and  Pratt  1979,  p.  Ill),  then  the  Gaussian  curvature  ne  is  the  product  /cj/c-y  and  the  mean 
curvature  Km  is  the  sum  ?q  +  k?.  For  a  surface  /(*,  y), 

,  (f\  _  [*jJuv  ~  fly 
8</J“(1+/2,  +  /2),‘ 


I 


h 


1  Which  is  why  Brady,  Grimson,  and  IitngridgWl9RC)  used  the  small  angle  approximation  f\x 
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Since  the  curvatures  can  be  negative  white  a  semi-noim  is  required  to  be  positive,  it  is  necessary  to 
investigate 


Grimson(1981)  observes  that  K*{af)  ^  |a| ic^(/)  because  of  the  denominator.  If  fx  and  fv  are  small,  the 
denominator  is  approximately  equal  to  one,  and  the  numerator  is  a  seminorm.  Note  that  it  is 


f**fw  /*y  (21) 

Grimson  shows  that  the  mean  curvature  /cm  is  also  not  a  semi-norm  for  exactly  the  same  reason.  The 
analogous  small  angle  approximation  is 


(/..  +  U2  =  (a/)j, 

the  square  Laplacian,  which  is  a  semi-norm.  We  find  it  convenient  to  denote  the  square  Laplacian  by  fj. 
Grimsont  1981)  chooses  the  quadratic  variation 


/L  +  2flv  +  fi 


w9 


on  the  grounds  that  its  null  space,  consisting  of  all  linear  functions,  is  smaller  thar  die  null  space  of  the  square 
Laplacian.  If  we  denote  the  quadratic  variation  by  Fq ,  we  sec  that  the  approximation  to  the  Gaussian  curvature 
given  in  equation  (21)  is 

How  shall  we  choose  a  performance  index  for  surface  interpolation,  given  that  it  has  to  satisfy  the  condi¬ 
tions  of  the  Corollary?  We  have  exhibited  three  candidates,  are  there  more?  Notice  first  that  each  of  the 
semi-norms  given  above  arc  quadratic  forms  in  /,„  /IV,  and  /„„.  It  is  easy  to  show  that  any  quadratic  form 
satisfies  the  semi-norm  and  parallelogram  conditions,  and  so  there  is  an  infinite  set  of  plausible  semi-norms  to 
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use  to  find  die  "most  conservative"  interpolated  surface.  We  need  an  extra  condition,  and  the  one  we  choose 
is  rotational  symmetry,  since  we  suppose  that  surface  interpolation  is  an  isotropic  process.  Proposition  6  of 
section  3  shows  that  the  rotationally  symmetric  quadratic  forms  in  /„,  /IW,  and  fvv  form  a  vector  space  that 
has  the  square  Laplacian  and  the  quadratic  variation  as  a  basis.  The  choice  of  which  performance  index  to  use 
is  thus  effectively  reduced  to  the  square  Laplacian,  the  quadratic  variation,  and  linear  combinations  of  them. 
How  shall  we  choose  between  those  two?  In  the  light  of  our  earlier  discussion,  two  criteria  suggest  themselves: 
the  Euler  equations  and  the  natural  boundary  conditions. 

Proposition  7.  All  rotationally  symmetric  quadratic  forms  lead  to  an  identical  Euler  equation 

Aa(/)  =  0. 

Proof.  We  exploit  the  fact  that  the  square  Laplacian  and  the  quadratic  variation  are  a  basis  of  the 
rotationally  symmetric  quadratic  forms. 

aSquare  Laplacian :  The  performance  index  is 


Ft  —  {fxz  -f/w)2- 


By  equation  (18)  the  Euler  equation  is 


+/„«  1}  +  JjsWA. + Av)}  “  ». 

that  is 

(A/)2  -  0, 


as  required. 

b .Quadratic  variation :  The  Kulcr  equation  is 
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2/ixxi  "f"  4ft yxy  "I"  2/j 


/VVW 


0, 


that  is 


(A/)2  =  0, 


provided  that  /  is  continuous  of  fourth  order. 

c.l. inear  combinations  offi  andFq. :  Linear  combinations  clearly  give  rise  to  the  identical  Euler  equation | 
The  gist  of  Proposition  7  is  that  there  is  no  difference  between  Fq  and  Fj  in  the  interior  away  from  the 
boundary  conditions.  We  can  see  the  result  of  Proposition  7  in  an  alternative  interesting  way.  Recall  that 


Fi—Fq  —  2{fxxfyy  —  fly), 

is  the  semi-norm  approximation  to  the  Gaussian  curvature  (equation  21).  The  latter  expression  is  an  instance 
of  a  divergence  expression ,  and  Courant  and  Hilbert(1953,  p.  1%)  note  "If  the  difference  between  the  in¬ 
tegrands  of  (wo  variational  problems  is  a  divergence  expression,  then  the  Euler  equations  and  therefore  the 
families  of  extremals  are  identical  for  the  two  variational  problems." 

Since  F,  and  F|  have  identical  Euler  equations,  we  analyze  their  natural  boundary  conditions  in  order 
to  choose  between  them.  We  could  approach  this  problem  directly;  but  a  more  revealing  route  is  available. 
Courant  and  Hilbcrt(1953,  p250)  consider  the  statics  of  a  thin  plate.  In  particular  they  determine  the  shape  it 
assumes  for  a  given  force  p(s)  along  its  boundary  T  and  bending  moments  m(a)  normal  to  its  boundary. 

Courant  and  Hilbert  note  that  the  energy  stored  in  the  plate  is  the  integral  of  a  quadratic  form  in  the 
principal  curvatures  «i  and  of  the  surface,  a  result  which  can  be  derived  from  noting  that  the  clastic  energy 
stored  in  a  thin  strip  (corresponding  to  any  normal  section)  is  proportional  to  the  square  curvature.  It  follows 
that  the  stored  energy  is  locally 
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6i  =  q(«j  -f  Kj)  -j-  2/3«i*i 

—  q(«i  4-  itaf  -f  2(/?  —  a>ci*i 
=  CMC*,  +  2(0  —  a)Kg 

**“'  aF(  2(0 

=  0F,  +  (a-0)F, 

=  o(MF,-}-(l-Mr,), 

a 

where  M  =  It  follows  that  the  energy  stored  in  the  thin  plate  is  a  convex  linear  combination  of  the 
square  Laplacian  and  the  quadratic  variation,  which  formally  establishes  its  connection  to  the  visual  percep¬ 
tual  problem  studied  here.  Observe  that  setting  the  weight/*  =  1  gives  the  square  Laplacian,  while  setting  it 
equal  t'-  nro  gives  the  quadratic  variation.  Note  also  that  this  expression  for  the  stored  energy  makes  use  of 
the  sm*'!  angle  approximation  to  the  curvature  used  in  equation  21. 

cond  source  of  stored  energy  derives  from  the  boundary  conditions  that  are  represented  as  a  function 
p(s)  a.  ,  the  boundary  T  of  the  plate  and  a  bending  moment  m(a)  applied  normal  to  the  plate.  Courant  and 
Hilbcrt(1953,  p.  251)  show  that  the  natural  boundary  conditions  associated  with  the  plate  are 

-A/-f  (1  -  n){fxxx]  +  2fxvxty,  +  fvvyl)  —  p{«) 

0  3 

—  A/  +  (1  —  rf^ifzxX'Xn  +  /xy(*,yn  +  Xn\J,)  +  /yyV.Vn)  —  ™(*)> 

that  is 


-A/  +  (1  -  v){\x.V']H\xtv,}T)  =  p(a) 

^  A/  -f  (1  -  ti)^(lxnyn]  H  \xty,]T)  =  m(a), 

where  H  is  the  Hessian  matrix 

fxx  fry 
fry  fyy 

Gladwell  and  Wait(  1979)  quote  version  of  this  result  due  to  Agmon(1965),  that  the  biharmonic  operator, 
which  we  showed  was  the  natural  boundary  condition  for  the  surface  interpolation  problem,  has  Dirichlct 
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forms  that  arc  linear  combinations  of  the  square  Laplacian  and  the  quadratic  variation.  As  an  example  of  the 
constraint,  consider  a  straight  line  contour  aligned  with  the  x-axis.  Then  [x,  yf]  =  [1 0]  and  [*„  y„]  =  [0 1]. 
The  natural  boundary  conditions  reduce  to 


fyy  4*  - P{*) 

fvvv  4“  (2  m)/w**  5=8 

The  constraint  is  tightest  when  n  is  not  equal  to  one.  A  similar  result  can  be  obtained  for  a  straight  line 
contour  inclined  at  an  angle  a  to  the  x-axis.  The  first  of  the  natural  boundary  conditions  is 


/n(sin2  a  4-  M  cos2  a)  4*  /vi/(«>s2  a  4-  M  sin2  a)  4~  (1  ~  m)  "in  2a/*,,. 

If  fi  —  1.  there  is  no  constraint  from  the  cross  derivative.  If  n  is  not  equal  to  1,  at  most  one  of  the 
terms  can  be  zero.  Wc  conclude  that  interpolation  problems  in  which  the  small  angle  approximations  used 
throughout  our  analysis  hold  it  is  preferable  to  choose  n  not  equal  to  one,  that  is  to  say  to  not  use  the  square 
1  .aplacian  as  a  performance  index.  The  quadratic  variation  is  an  obvious  choice,  but  so  are  linear  combinations 
of  the  square  1  .aplacian  and  the  quadratic  variation  for  which  n  is  not  equal  to  one.  Grimson(1981)  chooses 
the  quadratic  variation  since  its  null  space  is  smaller  than  that  of  the  square  laplacian.  litis  is  a  precise  way 
of  saying  that  it  imposes  a  tighter  constraint.  For  example,  the  function  /(x,  y)  —  xy  is  in  the  null  space  of 
the  square  Laplacian  but  not  in  the  null  space  of  the  quadratic  variation.  Since  the  quadratic  variation  has 
the  smallest  null  space  among  the  linear  combinations  of  the  square  laplacian  and  quadratic  variation,  this 
is  an  additional  reason  for  choosing  it.  We  "'ould  further  expect  that  any  differences  between  the  quadratic 
variation  and  the  square  Laplacian  would  snow  up  near  the  given  boundary  data  but  not  in  the  interior,  far 
removed  from  the  boundary.  Ihis  is  what  Grimson(1981)  finds  in  a  set  of  examples  that  compare  surfaces 
interpolated  using  the  quadialic  variation  and  the  square  Laplacian. 
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